Predicting catastrophes in nonlinear dynamical systems by compressive sensing 
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An extremely challenging problem of significant interest is to predict catastrophes in advance of their oc- 
currences. We present a general approach to predicting catastrophes in nonlinear dynamical systems under the 
assumption that the system equations are completely unknown and only time series reflecting the evolution of the 
dynamical variables of the system are available. Our idea is to expand the vector field or map of the underlying 
system into a suitable function series and then to use the compressive-sensing technique to accurately estimate 
the various terms in the expansion. Examples using paradigmatic chaotic systems are provided to demonstrate 
our idea. 
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PACS numbers; 05.45.-a 

It has been recognized that nonlinear dynamics are ubiqui- 
tous in many natural and engineering systems. A nonlinear 
system, in its parameter space, can often exhibit catastrophic 
bifurcations that ruin the desirable or "normal" state of oper- 
ation. Consider, for example, the phenomenon of crisis ^ 
where, as a system parameter is changed, a chaotic attrac- 
tor collides with its own basin boundary and is suddenly de- 
stroyed. After the crisis, the state of the system is completely 
different from that on the attractor before the crisis. Suppose 
that, for a nonlinear dynamical system, the state before the 
crisis is normal and desirable, and the state after the crisis is 
undesirable or destructive. The crisis can thus be regarded as 
a catastrophe that one strives to avoid at all cost. Catastrophic 
events, of course, can occur in different forms in all kinds 
of natural and man-made systems. A question of paramount 
importance is how to predict catastrophes in advance of their 
possible occurrences. This is especially challenging when the 
equations of the underlying dynamical system are unknown 
and one must then rely on measured time series or data to pre- 
dict any potential catastrophe. 

In this paper, we articulate a strategy to address the prob- 
lem of predicting catastrophes in nonlinear dynamical sys- 
tems. We assume that an accurate model of the system is not 
available, i.e., the system equations are unknown, but the time 
evolutions of the key variables of the system can be accessed 
through monitoring or measurements. Our method consists 
of three steps: (i) predicting the dynamical system based on 
time series, (ii) identifying the parameters of the system, and 
(iii) performing bifurcation analysis using the predicted sys- 
tem equations to locate potential catastrophic events in the pa- 
rameter space so as to determine the likelihood of system's 
drifting into a catastrophe regime. In particular, if the system 
operates at a parameter setting close to such a critical bifur- 
cation, catastrophe is imminent as a small parameter change 
or a random perturbation can push the system beyond the bi- 
furcation point. To be concrete, in this paper we regard crises 
as catastrophes. Once a complete set of system equations has 



been predicted and the parameters have been identified, one 
needs to examine the available parameter space. In general, to 
explore the multi-parameter space of a dynamical system can 
be extremely challenging, which can often lead to the discov- 
ery of new phenomena in dynamics. An early example in this 
area of research is the work by Stewart et al. |j2J, which in- 
vestigated the phenomenon of double crises in two-parameter 
dynamical systems. More recent efforts include the investi- 
gation of hierarchical structures in the parameter space dH. 
The present focus of our work, however, is on predicting the 
dynamical systems based on compressive sensing. 

The problem of predicting dynamical systems based on 
time series has been outstanding in nonlinear dynamics be- 
cause, despite previous efforts ||4|lin using the standard delay- 
coordinate embedding method f5l to decode the topological 
properties of the underlying system, how to accurately infer 
the underlying nonlinear system equations remains largely an 
unsolved problem. In principle, a nonlinear system can be ap- 
proximated by a large collection of linear equations in differ- 
ent regions of the phase space, which can indeed be achieved 
by reconstructing the Jacobian matrices on a proper grid that 
covers the phase-space region of interest [6, 7]. However, the 
accuracy and robustness of the procedure are challenging is- 
sues, including the difficulty with the required computations. 
In order to be able to predict potential catastrophes, local re- 
construction of a large set of linearized dynamics is not suffi- 
cient but rather, an accurate prediction of the underlying non- 
linear equations themselves is needed. 

Our framework to fully reconstruct dynamical systems us- 
ing time series alone is based on the assumption that the dy- 
namics of many natural and man-made systems are deter- 
mined by functions that can be approximated by series expan- 
sions in a suitable base. The major task is then to estimate the 
coefficients in the series representation. In general, the num- 
ber of coefficients to be estimated can be large but many of 
them are zero (the sparsity condition). According to the con- 
ventional wisdom this would be a difficult problem as a large 
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amount of data is required and the computations involved can 
be extremely demanding. However, the recent par adig m of 
compressive sensing developed by Candes et al. Js - 12 1 pro- 
vides a viable solution to the problem, where the key idea is to 
reconstruct a sparse signal from small amount of observations 
lH-lJ], as measured by linear projections of the original sig- 
nal on a few predetermined vectors. Since the requirements 
for the observations can be considerably relaxed as compared 
with those associated with conventional signal reconstruction 
schemes, compressive sensing has received much recent atten- 
tion and it is becoming a powerful technique to obtain high- 
fidelity signal for applications where sufficient observations 
are not available. Here, we shall articulate a general method- 
ology to cast the problems of dynamical-system prediction 
into the framework of compressive sensing and we demon- 
strate the power of our method by carrying out bifurcation 
analyses on the predicted dynamical systems to locate poten- 
tial catastrophes using exemplary chaotic systems. 

Generally, the problem of compressive sensing can be de- 
scribed as the reconstruction of a sparse vector a G R" from 
linear measurements X about a in the form: X = G ■ a, where 
X S R^, G is aw X V matrix and most components of a are 
zero. The compressive sensing theory ensures that the num- 
ber of components of the unknown signal can be much larger 
than the number of required measurements for reconstruc- 
tion, i.e., V ^ w. Accurate reconstruction can be achieved 
by solving the following convex optimization problem ist]: 
min ||a||i subject to X = G-a, where ||a||i = X^iLi l^^l 
is the Li norm of a. Solutions to the convex optimization 
problem have been worked out recently lH - 13 1 . 

We first show that the inverse problem of predicting dy- 
namical systems can be cast in the framework of compres- 
sive sensing so that optimal solutions can be obtained even 
when the number of base coefficients to be estimated is large 
and/or the amount of available data is small. In the follow- 
ing, we present a typical example to illustrate our method. 
Assume that the dynamical system can generally be written 
as X = F(x), where x e i?'" represents the set of externally 
accessible dynamical variables and F is a smooth vector func- 
tion in i?™. The jth component of F(x) can be represented as 
a power series: 
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where Xk {k = 1, ■ ■ ■ ,m) is the fcth component of the dy- 
namical variable, and the scalar coefficient of each product 
term (0^);^ ... e i? is to be determined from measure- 
ments. Note that the terms in Eq. ([T]l are all possible products 
of different components with different powers, and there are 
(1 + n)™ terms in total. 

To better explain our method, without loss of generality, we 
focus on one dynamical variable of the system. (Procedures 
for other variables are similar) For example, to construct 
the measurement vector X and the matrix G for the case of 
m = 3 (dynamical variables x, y, and z) and n = 3, we have 
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FIG. 1 : For the Henon map, in (a) x dimension and (b) y dimension, 
distributions of the predicted values of ten power-series coefficients 
up to order 3: constant, y, y , y , x, xy, xy , x , x y and x . 



the following explicit dynamical equation for the first dynam- 
ical variable: [F(x)]i = {ai)oM,ox°y° z° + {ai)i^o,ox^y° z° + 
■ ■ ■ + {ai)3.3,3X^y^z^- We can denote the coefficients of 
[F(x)]i by ai = [(ai)o,o.o, (ai)i,o,o, (01)3.3,3]^- 
Assuming that measurements of x(i) at a set of 
time ti,t2, ■ ■ ■ ,tiu are available, we denote g(i) = 
[xitfyitrzitf, x{tfy[tfz{t)\ • • • , x(t)^y{tf z{tf], 
such that [F(x(t))]i = g(t) • ai. From the expression 
of [F(x)]i, we can choose the measurement vector as 
X = [i(ti), i(t2), • • • which can be calculated 

from time series. Finally, we obtain the following equation in 
the form X = G • ai : 
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To ensure the restricted isometry property flSl, we normal- 
ize G by dividing elements in each column by the L2 norm 
of that column: (G')y = {G)ij/L2{j) with LzO') = 

Ej=i[(G)jj]2, so that X = G' • a.[. After the normal- 
ization, a'^ = ai ■ L2 can be determined via some standard 
compressive-sensing algorithm lil3ll . As a result, the coeffi- 
cients ai are given by sl^/L2- To determine the set of power- 
series coefficients corresponding to a different dynamical vari- 
able, say y, we simply replace the measurement vector by 
X = [y{ti),y{t2), ■ ■ ■ ,y{tw)]'^ and use the same matrix G. 
This way all coefficients can be estimated. 

We now present a number of physically relevant examples 
to illustrate our strategy. The first example is the Henon map 
1 1411 . a classical model that has been used to address many 
fundamental issues in chaotic dynamics. The prediction of 
map equations resembles that of a vector field. The map is 
given by: (x„+i, ?/„+i) = (1 - ax^ + ?/„, 6a;„), where a and 
6 are parameters. For h ~ 0.3, the map exhibits periodic and 
chaotic attractors for a < Uc ~ 1.42625, where is the crit- 
ical parameter value for a boundary crisis jUl, above which 
almost all trajectories diverge. The crisis can then be regarded 
as a catastrophe in the system evolution. Assuming, e.g., that 
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FIG. 2: Bifurcation diagrams of tiie predicted Henon map. The 
predicted map equations are: a;„+i — 0.999999996105743 + 
1.000000008610316 x?/„ ~axl andy„+i = 0.29999999837 x . 
The number Um of measurements used for prediction is 8 and the to- 
tal number n„z + nz of terms to be predicted is 16. 
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FIG. 3: Bifurcation diagram of (a) the predicted Lorenz system given 
by i = 10.000548307148881 x y - 10.001897147696283 x x, 
y = x(a- 1.000933186801829 x z) - 1.000092963203845 x y, 
i = 0.999893761636553 xa;i/ -2. 666477325955504 xz and of (b) 
the predicted Rossler system given by i = -0.999959701293536 x 
y - 0.999978902248041 x z, y = 1.000004981649221 x 
X + 0.200005996113158 y. y, i = 0.199997011085648 + 
0.999999156496251 x z{x - a). In both cases, n„ = 18 and 
+ nz = 35. 



the "normal" operation of the system corresponds to a chaotic 
attractor, we choose a = 1.4. Now suppose that the system 
operates at this parameter value and the system equations are 
completely unknown but the time series ({a;}„, {?/}«) can be 
obtained in real time. The goal is to assess, based on the time 
series only, how "close" the system is to a potential catastro- 
phe. (If measurements of only one dynamical variable can be 
obtained, one has to resort to the delay-coordinate embedding 
method |5].) For illustrative purpose, we assume power-series 
expansions up to order 3 in the map equations. Figure [T]shows 
the distributions of the estimated power-series coefficients, 
where we observe extremely narrow peaks about zero, indicat- 
ing that a large number of the coefficients are effectively zero, 
which correspond to nonexistent terms in the map equations. 
Coefficients that are not included in the zero peak coiTespond 
then to existent terms and they determine the predicted map 
equations. Figure |2] shows the bifurcation diagram from the 
predicted Henon map, which is consistent with the original 
diagram impressively well. In particular, the predicted sys- 
tem gives the value of the critical bifurcation point to within 
lO^'^, where a boundary crisis occur. Note that, to predict 
correctly the map equations, the number of required data is 
extremely low, not seen before in any method of dynamical- 
system reconstruction. Similar results have been obtained for 
the chaotic Lorenz [15J and Rossler oscillators, as shown 
by the predicted bifurcation diagrams in Figs. |3la) and|3lb), 
respectively. These agree with the original bifurcation dia- 
grams extremely well, so that all possible critical bifurcation 
points can be predicted accurately based on time series only. 

To quantify the performance of our method with respect to 
the amount of required data, we investigate the prediction er- 
rors which are defined separately for nonzero (existing) and 
zero terms in the dynamical equations. The relative eiTor of a 
nonzero term is defined as the ratio to the tiTie value of the ab- 
solute difference between the predicted and true values. The 
average over the errors of all terms in a component is the pre- 
diction error En^ of nonzero terms for the component. In 



contrast, the absolute error Ez is used for zero terms. Fig- 
ures Ufa) and Ub) show Enz as a function of the ratio of the 
number n,„ of measurements to the total number n„z + 
of terms to be predicted, for the standard map (vi\ and the 
Lorenz system, respectively. Note that, for the standard map, 
it is necessary to explore alternative bases of expansion so that 
the sparsity condition can be satisfied. Our strategy is that, as- 
suming a rough idea about the basic physics of the underlying 
dynamical system is available, we can choose a base that is 
compatible with the knowledge. In the case of the standard 
map, we thus choose the base which includes the trigonomet- 
ric functions. We obtain that, when the number n„i of mea- 
surements exceeds a threshold nt, Enz becomes effectively 
zero. Without loss of generality, we define rit by using the 
small threshold value E^z = lO^'' so that rit is the mini- 
mum number of required measurements for an accurate pre- 
diction. In Figs. Ua) andUb), we observe that rit is much 
less than n„2 + nz if Unz, the number of nonzero terms, is 
small. The performance of our method can thus be quantified 
by the threshold with respect to the numbers of measurements 
and terms to be predicted. As shown in Figs. lUc) andHtd) 
for the standard map and the Lorenz system, respectively, as 
the nonzero terms become sparser among all terms to be pre- 
dicted (characterized by a decrease in Unzjinnz + '^z) when 
^nz + f^z is increased), the ratio of the threshold to the 
total number of terms n^z + ^z becomes smaller. These re- 
sults demonstrate the advantage of our compressive-sensing 
based method to predict dynamical systems, i.e., high accu- 
racy and extremely low required measurements. In general, to 
predict the nonUnear dynamical system as accurately as pos- 
sible, many reasonable terms should be assumed in the expan- 
sions, insofar as the percentage of nonzero terms is small so 
that the sparsity condition of compressive sensing is satisfied. 

In addition, we examine the resistance of the method to 
measurement errors by inserting noise into time series. The 
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FIG. 4: Prediction errors Enz in dynamical equations as a function 
of the ratio of thie number Um of measurements to tlie total number 
n-nz + riz of terms to be predicted for (a) the standard map and (b) 
the Lorenz system. The ratio of the threshold rit to Unz + nz for 
different equations as a function of the ratio nnz/ijinz + riz) for (c) 
the standard map and (d) the Lorenz system. In (a) and (b), Unz + 
Uz is 20 and 35, respectively. The error bars represent the standard 
deviations obtained from 30 independent realizations. In (c) and (d), 
nnz + nz can be adjusted by the order of power series. In (c),the 
data points ranges from order 3 to order II, and in (d) from order 2 
to order 7. We find that Enz and Ez exhibit the same threshold. 



not work. A possible solution is to employ the Bayesian in- 
ference to determine the system equations. In general the 
computational challenge associated with the approach can be 
formidable, but the power-series or more general expansion 
based compressive-sensing method developed in this paper 
may present an effective strategy to overcome the difficulty. 

In summary, we have articulated a general approach to pre- 
dicting catastrophes in nonlinear dynamical systems. Our idea 
is to approximate the equations of the underlying system by 
series expansion and then to formulate the problem of esti- 
mating the various terms in the expansions using compressive 
sensing. The merit of our approach is then that, due to the 
nature of the compressive-sensing method, a large number of 
terms can be accurately estimated even with short available 
time series, enabling potential implementation in real times. 
We have presented a number of examples from chaotic dy- 
namics to demonstrate the effectiveness of our method. Pre- 
dicting catastrophe is a problem of uttermost importance in 
science and engineering and of extremely broad interest as 
well, and we hope our work will stimulate further efforts in 
this challenging area. 

This work was supported by AFOSR under Grants No. 
FA9550- 10- 1-0083 and FA9550-09- 1-0260. 
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FIG. 5: Prediction errors Enz as a function of noise amplitude for (a) 
the Henon map and (b) the standard map. Uniform noise is added to 
the time series. The values of rim and n„z + Uz for (a) are 8 and 16, 
respectively, and for (b) are 10 and 20, respectively. The prediction 
errors in the zero terms show similar behaviors. 



prediction errors as a function of noise amplitude are shown 
in Figs. 13 a) and|5jb) for the Henon map and the standard 
map, respectively. The results demonstrate that our method 
is robust against noise, due to the optimization nature of the 
compressive-sensing method. 

There are also situations where the system is high- 
dimensional or stochastic, for which the current method may 
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